Corona_Analysis

The 2019-2020 Coronavirus Pandemic Analysis

Contact: Smith Research

BACKGROUND & APPROACH

I wanted to track and trend the coronavirus outbreak on my own curiosity. There are some interesting questions that may fall out of this, as it is a very historic moment, including scientifically and analytically (we have a large amount of data being shared across the globe, analyzed in real-time). The world has come to a halt because of it.
This analysis attempts to answer the following questions (more to come):

1. What does the trend of the pandemic look like to date?
   
2. What are future case predictions based on historical model?
3. What interesting quirks or patterns emerge?

ASSUMPTIONS & LIMITATIONS: * This data is limited by the source. I realized early on that depending on source there were conflicting # of cases. Originally I was using JHU data… but this was always ‘ahead’ of the Our World In Data. I noticed that JHU’s website was buggy- you clicked on the U.S. stats but it didn’t reflect the U.S.. So I changed data sources to be more consistent with what is presented in the media (and Our World In Data has more extensive plots I can compare my own to). An interesting aside might be why the discrepancy? Was I missing something?
* Defintiions are important as is the idea that multiple varibales accumulate in things like total cases (more testing for example).

SOURCE RAW DATA: * https://ourworldindata.org/coronavirus
* https://github.com/CSSEGISandData/COVID-19/
*

INPUT DATA LOCATION: github (https://github.com/sbs87/coronavirus/tree/master/data)

OUTPUT DATA LOCATIOn: github (https://github.com/sbs87/coronavirus/tree/master/results)

TIMESTAMP

Start: ##—— Sun May 24 10:45:56 2020 ——##

PRE-ANALYSIS

The following sections are outside the scope of the ‘analysis’ but are still needed to prepare everything

UPSTREAM PROCESSING/ANALYSIS

1. Google Mobility Scraping, script available at get_google_mobility.py

# Mobility data has to be extracted from Google PDF reports using a web scraping script (python , written by Peter Simone, https://github.com/petersim1/MIT_COVID19)

# See get_google_mobility.py for local script 

python3 get_google_mobility.py
# writes csv file of mobility data as "mobility.csv"

SET UP ENVIORNMENT

Load libraries and set global variables

# timestamp start
timestamp()
## ##------ Sun May 24 10:45:56 2020 ------##

# clear previous enviornment
rm(list = ls())

##------------------------------------------
## LIBRARIES
##------------------------------------------
library(plyr)
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.0     ✓ purrr   0.3.3
## ✓ tibble  3.0.0     ✓ dplyr   0.8.5
## ✓ tidyr   1.0.2     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::arrange()   masks plyr::arrange()
## x purrr::compact()   masks plyr::compact()
## x dplyr::count()     masks plyr::count()
## x dplyr::failwith()  masks plyr::failwith()
## x dplyr::filter()    masks stats::filter()
## x dplyr::id()        masks plyr::id()
## x dplyr::lag()       masks stats::lag()
## x dplyr::mutate()    masks plyr::mutate()
## x dplyr::rename()    masks plyr::rename()
## x dplyr::summarise() masks plyr::summarise()
## x dplyr::summarize() masks plyr::summarize()
library(ggplot2)
library(reshape2)
## 
## Attaching package: 'reshape2'
## The following object is masked from 'package:tidyr':
## 
##     smiths
library(plot.utils)
library(utils)
library(knitr)

##------------------------------------------

##------------------------------------------
# GLOBAL VARIABLES
##------------------------------------------
user_name <- Sys.info()["user"]
working_dir <- paste0("/Users/", user_name, "/Projects/coronavirus/")  # don't forget trailing /
results_dir <- paste0(working_dir, "results/")  # assumes diretory exists
results_dir_custom <- paste0(results_dir, "custom/")  # assumes diretory exists


Corona_Cases.source_url <- "https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv"
Corona_Cases.US.source_url <- "https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_US.csv"
Corona_Deaths.US.source_url <- "https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_US.csv"
Corona_Deaths.source_url <- "https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_global.csv"

Corona_Cases.fn <- paste0(working_dir, "data/", basename(Corona_Cases.source_url))
Corona_Cases.US.fn <- paste0(working_dir, "data/", basename(Corona_Cases.US.source_url))
Corona_Deaths.fn <- paste0(working_dir, "data/", basename(Corona_Deaths.source_url))
Corona_Deaths.US.fn <- paste0(working_dir, "data/", basename(Corona_Deaths.US.source_url))
default_theme <- theme_bw() + theme(text = element_text(size = 14))  # fix this
##------------------------------------------

FUNCTIONS

List of functions

function_name	description
prediction_model	outputs case estumate for given log-linear moder parameters slope and intercept
make_long	converts input data to long format (specialized cases)
name_overlaps	outputs the column names intersection and set diffs of two data frame
find_linear_index	finds the first date at which linearaity occurs

##------------------------------------------
## FUNCTION: prediction_model
##------------------------------------------
## --- //// ----
# Takes days vs log10 (case) linear model parameters and a set of days since 100 cases and outputs a dataframe with total number of predicted cases for those days
## --- //// ----
prediction_model<-function(m=1,b=0,days=1){
  total_cases<-m*days+b
  total_cases.log<-log(total_cases,10)
  prediction<-data.frame(days=days,Total_confirmed_cases_perstate=total_cases)
  return(prediction)
}
##------------------------------------------

##------------------------------------------
## FUNCTION: make_long
##------------------------------------------
## --- //// ----
# Takes wide-format case data and converts into long format, using date and total cases as variable/values. Also enforces standardization/assumes data struture naming by using fixed variable name, value name, id.vars, 
## --- //// ----
make_long<-function(data_in,variable.name = "Date",
                   value.name = "Total_confirmed_cases",
                   id.vars=c("case_type","Province.State","Country.Region","Lat","Long","City","Population")){

long_data<-melt(data_in,
                id.vars = id.vars,
                variable.name=variable.name,
                value.name=value.name)
return(long_data)

}
##------------------------------------------

## THIS WILL BE IN UTILS AT SOME POINT
name_overlaps<-function(df1,df2){
i<-intersect(names(df1),
names(df2))
sd1<-setdiff(names(df1),
names(df2))
sd2<-setdiff(names(df2),names(df1))
cat("intersection:\n",paste(i,"\n"))
cat("in df1 but not df2:\n",paste(sd1,"\n"))
cat("in df2 but not df1:\n",paste(sd2,"\n"))
return(list("int"=i,"sd_1_2"=sd1,"sd_2_1"=sd2))
}

##------------------------------------------

##------------------------------------------
## FUNCTION: find_linear_index
##------------------------------------------
## --- //// ----
# Find date at which total case data is linear (for a given data frame) 
## --- //// ----

find_linear_index<-function(tmp,running_avg=5){
  tmp$Total_confirmed_cases_perstate.log<-log(tmp$Total_confirmed_cases_perstate,2)
  derivative<-data.frame(matrix(nrow = nrow(tmp),ncol = 4))
  names(derivative)<-c("m.time","mm.time","cumsum","date")
  
  # First derivative
  for(t in 2:nrow(tmp)){
    slope.t<- tmp[t,"Total_confirmed_cases_perstate.log"]- tmp[t-1,"Total_confirmed_cases_perstate.log"]
    derivative[t,"m.time"]<-slope.t
    derivative[t,"date"]<-as.Date(tmp[t,"Date"])
  }
  
  # Second derivative
  for(t in 2:nrow(derivative)){
    slope.t<- derivative[t,"m.time"]- derivative[t-1,"m.time"]
    derivative[t,"mm.time"]<-slope.t
  }
  
  #Compute running sum of second derivative (window = 5). Choose point at which within 0.2
  for(t in running_avg:nrow(derivative)){
    slope.t<- sum(abs(derivative[t:(t-4),"mm.time"])<0.2,na.rm = T)
    derivative[t,"cumsum"]<-slope.t
  }
  
  #Find date -5 from the stablility point
  linear_begin<-min(derivative[!is.na(derivative$cumsum) & derivative$cumsum==running_avg,"date"])-running_avg
  
  return(linear_begin)
}

READ IN DATA

• total number of cases. current source: https://github.com/CSSEGISandData (precvious source https://ourworldindata.org/coronavirus)

# Q: do we want to archive previous versions? Maybe an auto git mv?

##------------------------------------------
## Download and read in latest data from github
##------------------------------------------
download.file(Corona_Cases.source_url, destfile = Corona_Cases.fn)
Corona_Totals.raw <- read.csv(Corona_Cases.fn, header = T, stringsAsFactors = F)

download.file(Corona_Cases.US.source_url, destfile = Corona_Cases.US.fn)
Corona_Totals.US.raw <- read.csv(Corona_Cases.US.fn, header = T, stringsAsFactors = F)

download.file(Corona_Deaths.source_url, destfile = Corona_Deaths.fn)
Corona_Deaths.raw <- read.csv(Corona_Deaths.fn, header = T, stringsAsFactors = F)

download.file(Corona_Deaths.US.source_url, destfile = Corona_Deaths.US.fn)
Corona_Deaths.US.raw <- read.csv(Corona_Deaths.US.fn, header = T, stringsAsFactors = F)

# latest date on all data:
paste("US deaths:", names(Corona_Deaths.US.raw)[ncol(Corona_Deaths.US.raw)])

## [1] "US deaths: X5.23.20"

paste("US total:", names(Corona_Totals.US.raw)[ncol(Corona_Totals.US.raw)])

## [1] "US total: X5.23.20"

paste("World deaths:", names(Corona_Deaths.raw)[ncol(Corona_Deaths.raw)])

## [1] "World deaths: X5.23.20"

paste("World total:", names(Corona_Totals.raw)[ncol(Corona_Totals.raw)])

## [1] "World total: X5.23.20"

PROCESS DATA

• Convert to long format
  
• Fix date formatting/convert to numeric date
  
• Log10 transform total # cases

##------------------------------------------
## Combine death and total data frames
##------------------------------------------
Corona_Totals.raw$case_type<-"total"
Corona_Totals.US.raw$case_type<-"total"
Corona_Deaths.raw$case_type<-"death"
Corona_Deaths.US.raw$case_type<-"death"

# for some reason, Population listed in US death file but not for other data... Weird. When combining, all datasets will have this column, but US deaths is the only useful one.  
Corona_Totals.US.raw$Population<-"NA" 
Corona_Totals.raw$Population<-"NA"
Corona_Deaths.raw$Population<-"NA"

Corona_Cases.raw<-rbind(Corona_Totals.raw,Corona_Deaths.raw)
Corona_Cases.US.raw<-rbind(Corona_Totals.US.raw,Corona_Deaths.US.raw)
#TODO: custom utils- setdiff, intersect names... option to output in merging too
##------------------------------------------
# prepare raw datasets for eventual combining
##------------------------------------------
Corona_Cases.raw$City<-"NA" # US-level data has Cities
Corona_Cases.US.raw$Country_Region<-"US_state" # To differentiate from World-level stats

Corona_Cases.US.raw<-plyr::rename(Corona_Cases.US.raw,c("Province_State"="Province.State",
                                                  "Country_Region"="Country.Region",
                                                  "Long_"="Long",
                                                  "Admin2"="City"))


##------------------------------------------
## Convert to long format
##------------------------------------------
#JHU has a gross file format. It's in wide format with each column is the date in MM/DD/YY. So read this in as raw data but trasnform it to be better suited for analysis
# Furthermore, the World and US level data is formatted differently, containing different columns, etc. Recitfy this and combine the world-level stats with U.S. level stats.

Corona_Cases.long<-rbind(make_long(select(Corona_Cases.US.raw,-c(UID,iso2,iso3,code3,FIPS,Combined_Key))),
make_long(Corona_Cases.raw))


##------------------------------------------
## Fix date formatting, convert to numeric date
##------------------------------------------
Corona_Cases.long$Date<-gsub(Corona_Cases.long$Date,pattern = "^X",replacement = "0") # leading 0 read in as X
Corona_Cases.long$Date<-gsub(Corona_Cases.long$Date,pattern = "20$",replacement = "2020") # ends in .20 and not 2020
Corona_Cases.long$Date<-as.Date(Corona_Cases.long$Date,format = "%m.%d.%y")
Corona_Cases.long$Date.numeric<-as.numeric(Corona_Cases.long$Date)

kable(table(select(Corona_Cases.long,c("Country.Region","case_type"))),caption = "Number of death and total case longitudinal datapoints per geographical region")

Number of death and total case longitudinal datapoints per geographical region

	death	total
Afghanistan	123	123
Albania	123	123
Algeria	123	123
Andorra	123	123
Angola	123	123
Antigua and Barbuda	123	123
Argentina	123	123
Armenia	123	123
Australia	984	984
Austria	123	123
Azerbaijan	123	123
Bahamas	123	123
Bahrain	123	123
Bangladesh	123	123
Barbados	123	123
Belarus	123	123
Belgium	123	123
Belize	123	123
Benin	123	123
Bhutan	123	123
Bolivia	123	123
Bosnia and Herzegovina	123	123
Botswana	123	123
Brazil	123	123
Brunei	123	123
Bulgaria	123	123
Burkina Faso	123	123
Burma	123	123
Burundi	123	123
Cabo Verde	123	123
Cambodia	123	123
Cameroon	123	123
Canada	1722	1722
Central African Republic	123	123
Chad	123	123
Chile	123	123
China	4059	4059
Colombia	123	123
Comoros	123	123
Congo (Brazzaville)	123	123
Congo (Kinshasa)	123	123
Costa Rica	123	123
Cote d’Ivoire	123	123
Croatia	123	123
Cuba	123	123
Cyprus	123	123
Czechia	123	123
Denmark	369	369
Diamond Princess	123	123
Djibouti	123	123
Dominica	123	123
Dominican Republic	123	123
Ecuador	123	123
Egypt	123	123
El Salvador	123	123
Equatorial Guinea	123	123
Eritrea	123	123
Estonia	123	123
Eswatini	123	123
Ethiopia	123	123
Fiji	123	123
Finland	123	123
France	1353	1353
Gabon	123	123
Gambia	123	123
Georgia	123	123
Germany	123	123
Ghana	123	123
Greece	123	123
Grenada	123	123
Guatemala	123	123
Guinea	123	123
Guinea-Bissau	123	123
Guyana	123	123
Haiti	123	123
Holy See	123	123
Honduras	123	123
Hungary	123	123
Iceland	123	123
India	123	123
Indonesia	123	123
Iran	123	123
Iraq	123	123
Ireland	123	123
Israel	123	123
Italy	123	123
Jamaica	123	123
Japan	123	123
Jordan	123	123
Kazakhstan	123	123
Kenya	123	123
Korea, South	123	123
Kosovo	123	123
Kuwait	123	123
Kyrgyzstan	123	123
Laos	123	123
Latvia	123	123
Lebanon	123	123
Lesotho	123	123
Liberia	123	123
Libya	123	123
Liechtenstein	123	123
Lithuania	123	123
Luxembourg	123	123
Madagascar	123	123
Malawi	123	123
Malaysia	123	123
Maldives	123	123
Mali	123	123
Malta	123	123
Mauritania	123	123
Mauritius	123	123
Mexico	123	123
Moldova	123	123
Monaco	123	123
Mongolia	123	123
Montenegro	123	123
Morocco	123	123
Mozambique	123	123
MS Zaandam	123	123
Namibia	123	123
Nepal	123	123
Netherlands	615	615
New Zealand	123	123
Nicaragua	123	123
Niger	123	123
Nigeria	123	123
North Macedonia	123	123
Norway	123	123
Oman	123	123
Pakistan	123	123
Panama	123	123
Papua New Guinea	123	123
Paraguay	123	123
Peru	123	123
Philippines	123	123
Poland	123	123
Portugal	123	123
Qatar	123	123
Romania	123	123
Russia	123	123
Rwanda	123	123
Saint Kitts and Nevis	123	123
Saint Lucia	123	123
Saint Vincent and the Grenadines	123	123
San Marino	123	123
Sao Tome and Principe	123	123
Saudi Arabia	123	123
Senegal	123	123
Serbia	123	123
Seychelles	123	123
Sierra Leone	123	123
Singapore	123	123
Slovakia	123	123
Slovenia	123	123
Somalia	123	123
South Africa	123	123
South Sudan	123	123
Spain	123	123
Sri Lanka	123	123
Sudan	123	123
Suriname	123	123
Sweden	123	123
Switzerland	123	123
Syria	123	123
Taiwan*	123	123
Tajikistan	123	123
Tanzania	123	123
Thailand	123	123
Timor-Leste	123	123
Togo	123	123
Trinidad and Tobago	123	123
Tunisia	123	123
Turkey	123	123
Uganda	123	123
Ukraine	123	123
United Arab Emirates	123	123
United Kingdom	1353	1353
Uruguay	123	123
US	123	123
US_state	401103	401103
Uzbekistan	123	123
Venezuela	123	123
Vietnam	123	123
West Bank and Gaza	123	123
Western Sahara	123	123
Yemen	123	123
Zambia	123	123
Zimbabwe	123	123

# Decouple population and lat/long data, refactor to make it more tidy
metadata_columns<-c("Lat","Long","Population")
metadata<-unique(select(filter(Corona_Cases.long,case_type=="death"),c("Country.Region","Province.State","City",all_of(metadata_columns))))
Corona_Cases.long<-select(Corona_Cases.long,-all_of(metadata_columns))

# Some counties are not summarized on the country level. collapse all but US
Corona_Cases.long<-rbind.fill(ddply(filter(Corona_Cases.long,!Country.Region=="US_state"),c("case_type","Country.Region","Date","Date.numeric"),summarise,Total_confirmed_cases=sum(Total_confirmed_cases)),filter(Corona_Cases.long,Country.Region=="US_state"))

# Put total case and deaths side-by-side (wide)
Corona_Cases<-spread(Corona_Cases.long,key = case_type,value = Total_confirmed_cases)

#Compute moratlity rate
Corona_Cases$mortality_rate<-Corona_Cases$death/Corona_Cases$total

#TMP
Corona_Cases<-plyr::rename(Corona_Cases,c("total"="Total_confirmed_cases","death"="Total_confirmed_deaths"))

##------------------------------------------
## log10 transform total # cases
##------------------------------------------
Corona_Cases$Total_confirmed_cases.log<-log(Corona_Cases$Total_confirmed_cases,10)
Corona_Cases$Total_confirmed_deaths.log<-log(Corona_Cases$Total_confirmed_deaths,10)
##------------------------------------------
       
##------------------------------------------
## Compute # of days since 100th for US data
##------------------------------------------

# Find day that 100th case was found for Country/Province. NOTE: Non US countries may have weird provinces. For example, Frane is summairzed at the country level but also had 3 providences. I've only ensured the U.S. case100 works... so the case100_date for U.S. is summarized both for the entire country (regardless of state) and on a per-state level. 
# TODO: consider city-level summary as well. This data may be sparse

Corona_Cases<-merge(Corona_Cases,ddply(filter(Corona_Cases,Total_confirmed_cases>100),c("Country.Region"),summarise,case100_date=min(Date.numeric)))
Corona_Cases$Days_since_100<-Corona_Cases$Date.numeric-Corona_Cases$case100_date

##------------------------------------------
## Add population and lat/long data (CURRENTLY US ONLY)
##------------------------------------------

kable(filter(metadata,(is.na(Country.Region) | is.na(Population) )) %>% select(c("Country.Region","Province.State","City")) %>% unique(),caption = "Regions for which either population or Country is NA")

Regions for which either population or Country is NA

Country.Region	Province.State	City

# Drop missing data 
metadata<-filter(metadata,!(is.na(Country.Region) | is.na(Population) ))
# Convert remaining pop to numeric
metadata$Population<-as.numeric(metadata$Population)

## Warning: NAs introduced by coercion

# Add metadata to cases
Corona_Cases<-merge(Corona_Cases,metadata,all.x = T)

##------------------------------------------
## Compute total and death cases relative to population 
##------------------------------------------

Corona_Cases$Total_confirmed_cases.per100<-100*Corona_Cases$Total_confirmed_cases/Corona_Cases$Population
Corona_Cases$Total_confirmed_deaths.per100<-100*Corona_Cases$Total_confirmed_deaths/Corona_Cases$Population


##------------------------------------------
## Filter df for US state-wide stats
##------------------------------------------

Corona_Cases.US_state<-filter(Corona_Cases,Country.Region=="US_state" & Total_confirmed_cases>0 ) 
kable(table(select(Corona_Cases.US_state,c("Province.State"))),caption = "Number of longitudinal datapoints (total/death) per state")

Number of longitudinal datapoints (total/death) per state

Var1	Freq
Alabama	4065
Alaska	747
Arizona	1031
Arkansas	4230
California	3933
Colorado	3706
Connecticut	610
Delaware	251
Diamond Princess	68
District of Columbia	69
Florida	4362
Georgia	9688
Grand Princess	69
Guam	69
Hawaii	357
Idaho	1925
Illinois	5502
Indiana	5641
Iowa	5085
Kansas	4173
Kentucky	6119
Louisiana	4085
Maine	1035
Maryland	1601
Massachusetts	1043
Michigan	4915
Minnesota	4588
Mississippi	5062
Missouri	5538
Montana	1730
Nebraska	3114
Nevada	765
New Hampshire	719
New Jersey	1553
New Mexico	1684
New York	3917
North Carolina	5936
North Dakota	1966
Northern Mariana Islands	54
Ohio	5301
Oklahoma	3976
Oregon	2042
Pennsylvania	4199
Puerto Rico	69
Rhode Island	406
South Carolina	2961
South Dakota	2509
Tennessee	5719
Texas	11862
Utah	954
Vermont	960
Virgin Islands	69
Virginia	7462
Washington	2683
West Virginia	2694
Wisconsin	3986
Wyoming	1231

Corona_Cases.US_state<-merge(Corona_Cases.US_state,ddply(filter(Corona_Cases.US_state,Total_confirmed_cases>100),c("Province.State"),summarise,case100_date_state=min(Date.numeric)))
Corona_Cases.US_state$Days_since_100_state<-Corona_Cases.US_state$Date.numeric-Corona_Cases.US_state$case100_date_state

ANALYSIS

Q1: What is the trend in cases, mortality across geopgraphical regions?

Plot # of cases vs time
* For each geographical set:
* comparative longitudinal case trend (absolute & log scale)
* comparative longitudinal mortality trend
* death vs total correlation

question	dataset	x	y	color	facet	pch	dimentions
comparative_longitudinal_case_trend	long	time	log_cases	geography	none (case type?)	case_type	[15, 50, 4] geography x (2 scale?) case type
comparative longitudinal case trend	long	time	cases	geography	case_type	?	[15, 50, 4] geography x (2+ scale) case type
comparative longitudinal mortality trend	wide	time	mortality rate	geography	none	none	[15, 50, 4] geography
death vs total correlation	wide	cases	deaths	geography	none	none	[15, 50, 4] geography

# total cases vs time
# death cases vs time
# mortality rate vs time
# death vs mortality


  # death vs mortality
  # total & death case vs time (same plot)

#<question> <x> <y> <colored> <facet> <dataset>
## trend in case/deaths over time, comapred across regions <time> <log cases> <geography*> <none> <.wide>
## trend in case/deaths over time, comapred across regions <time> <cases> <geography*> <case_type> <.long>
## trend in mortality rate over time, comapred across regions <time> <mortality rate> <geography*> <none>
## how are death/mortality related/correlated? <time> <log cases> <geography*> <none>
## how are death and case load correlated? <cases> <deaths>

# lm for each?? - > apply lm from each region starting from 100th case. m, b associated with each.
    # input: geographical regsion, logcase vs day (100th case)
    # output: m, b for each geographical region ID



#total/death on same plot-  diffeer by 2 logs, so when plotting log, use pch. when plotting absolute, need to use free scales
#when plotting death and case on same, melt. 

#CoronaCases - > filter sets (3)
  #world - choose countries with sufficent data

N<-ddply(filter(Corona_Cases,Total_confirmed_cases>100),c("Country.Region"),summarise,n=length(Country.Region))
ggplot(filter(N,n<100),aes(x=n))+
  geom_histogram()+
  default_theme+
  ggtitle("Distribution of number of days with at least 100 confirmed cases for each region")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

kable(arrange(N,-n),caption="Sorted number of days with at least 100 confirmed cases")

Sorted number of days with at least 100 confirmed cases

Country.Region	n
US_state	35562
China	123
Diamond Princess	104
Korea, South	94
Japan	93
Italy	91
Iran	88
Singapore	85
France	84
Germany	84
Spain	83
US	82
Switzerland	80
United Kingdom	80
Belgium	79
Netherlands	79
Norway	79
Sweden	79
Austria	77
Malaysia	76
Australia	75
Bahrain	75
Denmark	75
Canada	74
Qatar	74
Iceland	73
Brazil	72
Czechia	72
Finland	72
Greece	72
Iraq	72
Israel	72
Portugal	72
Slovenia	72
Egypt	71
Estonia	71
India	71
Ireland	71
Kuwait	71
Philippines	71
Poland	71
Romania	71
Saudi Arabia	71
Indonesia	70
Lebanon	70
San Marino	70
Thailand	70
Chile	69
Pakistan	69
Luxembourg	68
Peru	68
Russia	68
Ecuador	67
Mexico	67
Slovakia	67
South Africa	67
United Arab Emirates	67
Armenia	66
Colombia	66
Croatia	66
Panama	66
Serbia	66
Taiwan*	66
Turkey	66
Argentina	65
Bulgaria	65
Latvia	65
Uruguay	65
Algeria	64
Costa Rica	64
Dominican Republic	64
Hungary	64
Andorra	63
Bosnia and Herzegovina	63
Jordan	63
Lithuania	63
Morocco	63
New Zealand	63
North Macedonia	63
Vietnam	63
Albania	62
Cyprus	62
Malta	62
Moldova	62
Brunei	61
Burkina Faso	61
Sri Lanka	61
Tunisia	61
Ukraine	60
Azerbaijan	59
Ghana	59
Kazakhstan	59
Oman	59
Senegal	59
Venezuela	59
Afghanistan	58
Cote d’Ivoire	58
Cuba	57
Mauritius	57
Uzbekistan	57
Cambodia	56
Cameroon	56
Honduras	56
Nigeria	56
West Bank and Gaza	56
Belarus	55
Georgia	55
Bolivia	54
Kosovo	54
Kyrgyzstan	54
Montenegro	54
Congo (Kinshasa)	53
Kenya	52
Niger	51
Guinea	50
Rwanda	50
Trinidad and Tobago	50
Paraguay	49
Bangladesh	48
Djibouti	46
El Salvador	45
Guatemala	44
Madagascar	43
Mali	42
Congo (Brazzaville)	39
Jamaica	39
Gabon	37
Somalia	37
Tanzania	37
Ethiopia	36
Burma	35
Sudan	34
Liberia	33
Maldives	31
Equatorial Guinea	30
Cabo Verde	28
Sierra Leone	26
Guinea-Bissau	25
Togo	25
Zambia	24
Eswatini	23
Chad	22
Tajikistan	21
Haiti	19
Sao Tome and Principe	19
Benin	17
Nepal	17
Uganda	17
Central African Republic	16
South Sudan	16
Guyana	14
Mozambique	13
Yemen	9
Mongolia	8
Mauritania	5
Nicaragua	5

# Pick top 15 countries with data
max_colors<-12
# find way to fix this- China has diff provences. Plot doesnt look right...
sufficient_data<-arrange(filter(N,!Country.Region %in% c("US_state", "Diamond Princess")),-n)[1:max_colors,]
kable(sufficient_data,caption = paste0("Top ",max_colors," countries with sufficient data"))

Top 12 countries with sufficient data

Country.Region	n
China	123
Korea, South	94
Japan	93
Italy	91
Iran	88
Singapore	85
France	84
Germany	84
Spain	83
US	82
Switzerland	80
United Kingdom	80

Corona_Cases.world<-filter(Corona_Cases,Country.Region %in% c(sufficient_data$Country.Region))


  #us 
  #    - by state
Corona_Cases.US<-filter(Corona_Cases,Country.Region=="US" & Total_confirmed_cases>0)
# summarize 
#!City %in% c("Unassigned") 
  #    - specific cities
#mortality_rate!=Inf & mortality_rate<=1
Corona_Cases.UScity<-filter(Corona_Cases,Province.State %in% c("Pennsylvania","Maryland","New York","New Jersey") & City %in% c("Bucks","Baltimore City", "New York","Burlington","Cape May"))

measure_vars_long<-c("Total_confirmed_cases.log","Total_confirmed_cases","Total_confirmed_deaths","Total_confirmed_deaths.log")
melt_arg_list<-list(variable.name = "case_type",value.name = "cases",measure.vars = c("Total_confirmed_cases","Total_confirmed_deaths"))
melt_arg_list$data=NULL


melt_arg_list$data=select(Corona_Cases.world,-ends_with(match = "log"))
Corona_Cases.world.long<-do.call(melt,melt_arg_list)
melt_arg_list$data=select(Corona_Cases.UScity,-ends_with(match = "log"))
Corona_Cases.UScity.long<-do.call(melt,melt_arg_list)
melt_arg_list$data=select(Corona_Cases.US_state,-ends_with(match = "log"))
Corona_Cases.US_state.long<-do.call(melt,melt_arg_list)

Corona_Cases.world.long$cases.log<-log(Corona_Cases.world.long$cases,10)
Corona_Cases.US_state.long$cases.log<-log(Corona_Cases.US_state.long$cases,10)
Corona_Cases.UScity.long$cases.log<-log(Corona_Cases.UScity.long$cases,10)


# what is the current death and total case load for US? For world? For states?
#-absolute
#-log

# what is mortality rate (US, world)
#-absolute

#how is death and case correlated? (US, world)
#-absolute

#Corona_Cases.US<-filter(Corona_Cases,Country.Region=="US" & Total_confirmed_cases>0)
#Corona_Cases.US.case100<-filter(Corona_Cases.US, Days_since_100>=0)
# linear model parameters
#(model_fit<-lm(formula = Total_confirmed_cases.log~Days_since_100,data= Corona_Cases.US.case100 ))

#(slope<-model_fit$coefficients[2])
#(intercept<-model_fit$coefficients[1])

# Correlation coefficient
#cor(x = Corona_Cases.US.case100$Days_since_100,y = Corona_Cases.US.case100$Total_confirmed_cases.log)

##------------------------------------------
## Plot World Data
##------------------------------------------
# Timestamp for world
timestamp_plot.world<-paste("Most recent date for which data available:",max(Corona_Cases.world$Date))#timestamp(quiet = T,prefix = "Updated ",suffix = " (EST)")


# Base template for plots
baseplot.world<-ggplot(data=NULL,aes(x=Days_since_100,col=Country.Region))+
  default_theme+
  scale_color_brewer(type = "qualitative",palette = "Paired")+
  ggtitle(paste("Log10 cases over time,",timestamp_plot.world))+
  theme(legend.position = "bottom",plot.title = element_text(size=12))


##/////////////////////////
### Plot Longitudinal cases

(Corona_Cases.world.long.plot<-baseplot.world+
    geom_point(data=Corona_Cases.world.long,aes(y=cases))+
    geom_line(data=Corona_Cases.world.long,aes(y=cases))+
    facet_wrap(~case_type,scales = "free_y",ncol=1)+
    ggtitle(timestamp_plot.world)
    )

(Corona_Cases.world.loglong.plot<-baseplot.world+
    geom_point(data=Corona_Cases.world.long,aes(y=cases.log))+
    geom_line(data=Corona_Cases.world.long,aes(y=cases.log))+
    facet_wrap(~case_type,scales = "free_y",ncol=1)+
    ggtitle(timestamp_plot.world))

##/////////////////////////
### Plot Longitudinal mortality rate

(Corona_Cases.world.mortality.plot<-baseplot.world+
    geom_point(data=Corona_Cases.world,aes(y=mortality_rate))+
    geom_line(data=Corona_Cases.world,aes(y=mortality_rate))+
    ylim(c(0,0.3))+
    ggtitle(timestamp_plot.world))

## Warning: Removed 100 rows containing missing values (geom_point).

## Warning: Removed 100 row(s) containing missing values (geom_path).

##/////////////////////////
### Plot death vs total case correlation

(Corona_Cases.world.casecor.plot<-ggplot(Corona_Cases.world,aes(x=Total_confirmed_cases,y=Total_confirmed_deaths,col=Country.Region))+
  geom_point()+
  geom_line()+
  default_theme+
  scale_color_brewer(type = "qualitative",palette = "Paired")+
  ggtitle(paste("Log10 cases over time,",timestamp_plot.world))+
  theme(legend.position = "bottom",plot.title = element_text(size=12))+
    ggtitle(timestamp_plot.world))

### Write polots

write_plot(Corona_Cases.world.long.plot,wd = results_dir)

## [1] "/Users/stevensmith/Projects/coronavirus/results/Corona_Cases.world.long.plot.png"

write_plot(Corona_Cases.world.loglong.plot,wd = results_dir)

## [1] "/Users/stevensmith/Projects/coronavirus/results/Corona_Cases.world.loglong.plot.png"

write_plot(Corona_Cases.world.mortality.plot,wd = results_dir)

## Warning: Removed 100 rows containing missing values (geom_point).

## Warning: Removed 100 row(s) containing missing values (geom_path).

## [1] "/Users/stevensmith/Projects/coronavirus/results/Corona_Cases.world.mortality.plot.png"

write_plot(Corona_Cases.world.casecor.plot,wd = results_dir)

## [1] "/Users/stevensmith/Projects/coronavirus/results/Corona_Cases.world.casecor.plot.png"

##------------------------------------------
## Plot US State Data
##-----------------------------------------

baseplot.US<-ggplot(data=NULL,aes(x=Days_since_100_state,col=case_type))+
  default_theme+
  facet_wrap(~Province.State)+
  ggtitle(paste("Log10 cases over time,",timestamp_plot.world))

Corona_Cases.US_state.long.plot<-baseplot.US+geom_point(data=Corona_Cases.US_state.long,aes(y=cases.log))
##------------------------------------------
## Plot US City Data
##-----------------------------------------

Corona_Cases.US.plotdata<-filter(Corona_Cases.US_state,Province.State %in% c("Pennsylvania","Maryland","New York","New Jersey") &
                                   City %in% c("Bucks","Baltimore City", "New York","Burlington","Cape May") &
                                   Total_confirmed_cases>0) 
timestamp_plot<-paste("Most recent date for which data available:",max(Corona_Cases.US.plotdata$Date))#timestamp(quiet = T,prefix = "Updated ",suffix = " (EST)")

city_colors<-c("Bucks"='#beaed4',"Baltimore City"='#386cb0', "New York"='#7fc97f',"Burlington"='#fdc086',"Cape May"="#e78ac3")

##/////////////////////////
### Plot death vs total case correlation

(Corona_Cases.city.loglong.plot<-ggplot(melt(Corona_Cases.US.plotdata,measure.vars = c("Total_confirmed_cases.log","Total_confirmed_deaths.log"),variable.name = "case_type",value.name = "cases"),aes(x=Date,y=cases,col=City,pch=case_type))+
  geom_point(size=4)+
    geom_line()+
  default_theme+
  #facet_wrap(~case_type)+
    ggtitle(paste("Log10 total and death cases over time,",timestamp_plot))+
theme(legend.position = "bottom",plot.title = element_text(size=12),axis.text.x = element_text(angle=45,hjust=1))+
    scale_color_manual(values = city_colors)+
  scale_x_date(date_breaks="1 week",date_minor_breaks="1 day"))

(Corona_Cases.city.long.plot<-ggplot(filter(Corona_Cases.US.plotdata,Province.State !="New York"),aes(x=Date,y=Total_confirmed_cases,col=City))+
  geom_point(size=4)+
  geom_line()+
  default_theme+
  facet_grid(~Province.State,scales = "free_y")+
  ggtitle(paste("MD, PA, NJ total cases over time,",timestamp_plot))+
  theme(legend.position = "bottom",plot.title = element_text(size=12),axis.text.x = element_text(angle=45,hjust=1))
+
  scale_color_manual(values = city_colors)+
  scale_x_date(date_breaks="1 week",date_minor_breaks="1 day"))

(Corona_Cases.city.mortality.plot<-ggplot(Corona_Cases.US.plotdata,aes(x=Date,y=mortality_rate,col=City))+
  geom_point(size=3)+
  geom_line(size=2)+
  default_theme+
  ggtitle(paste("Mortality rate (deaths/total) over time,",timestamp_plot))+
  theme(legend.position = "bottom",plot.title = element_text(size=12),axis.text.x = element_text(angle=45,hjust=1))+
  scale_color_manual(values = city_colors)+
  scale_x_date(date_breaks="1 week",date_minor_breaks="1 day"))

(Corona_Cases.city.casecor.plot<-ggplot(filter(Corona_Cases.US.plotdata,Province.State !="New York"),aes(y=Total_confirmed_deaths,x=Total_confirmed_cases,col=City))+
  geom_point(size=3)+
  geom_line(size=2)+
  default_theme+
  ggtitle(paste("Correlation of death vs total cases,",timestamp_plot))+
  theme(legend.position = "bottom",plot.title = element_text(size=12))+
  scale_color_manual(values = city_colors))

(Corona_Cases.city.long.normalized.plot<-ggplot(filter(Corona_Cases.US.plotdata,Province.State !="New York"),aes(x=Date,y=Total_confirmed_cases.per100,col=City))+
  geom_point(size=4)+
  geom_line()+
  default_theme+
  facet_grid(~Province.State)+
  ggtitle(paste("MD, PA, NJ total cases over time per 100 people,",timestamp_plot))+
  theme(legend.position = "bottom",plot.title = element_text(size=12),axis.text.x = element_text(angle=45,hjust=1))+
  scale_color_manual(values = city_colors)  +
  scale_x_date(date_breaks="1 week",date_minor_breaks="1 day"))

write_plot(Corona_Cases.city.long.plot,wd = results_dir_custom)

## [1] "/Users/stevensmith/Projects/coronavirus/results/custom/Corona_Cases.city.long.plot.png"

write_plot(Corona_Cases.city.loglong.plot,wd = results_dir_custom)

## [1] "/Users/stevensmith/Projects/coronavirus/results/custom/Corona_Cases.city.loglong.plot.png"

write_plot(Corona_Cases.city.mortality.plot,wd = results_dir_custom)

## [1] "/Users/stevensmith/Projects/coronavirus/results/custom/Corona_Cases.city.mortality.plot.png"

write_plot(Corona_Cases.city.casecor.plot,wd = results_dir_custom)

## [1] "/Users/stevensmith/Projects/coronavirus/results/custom/Corona_Cases.city.casecor.plot.png"

write_plot(Corona_Cases.city.long.normalized.plot,wd = results_dir_custom)

## [1] "/Users/stevensmith/Projects/coronavirus/results/custom/Corona_Cases.city.long.normalized.plot.png"

Q1b what is the model

Fit the cases to a linear model 1. Find time at which the case vs date becomes linear in each plot
2. Fit linear model for each city

# What is the predict # of cases for the next few days?
# How is the model performing historically?

Corona_Cases.US_state.summary<-ddply(Corona_Cases.US_state,
                                     c("Province.State","Date"),
                                     summarise,
                                     Total_confirmed_cases_perstate=sum(Total_confirmed_cases)) %>% 
    filter(Total_confirmed_cases_perstate>100)

# Compute the states with the most cases (for coloring and for linear model)
top_states_totals<-head(ddply(Corona_Cases.US_state.summary,c("Province.State"),summarise, Total_confirmed_cases_perstate.max=max(Total_confirmed_cases_perstate)) %>% arrange(-Total_confirmed_cases_perstate.max),n=max_colors)

kable(top_states_totals,caption = "Top 12 States, total count ")
top_states<-top_states_totals$Province.State

# Manually fix states so that Maryland is switched out for New York
top_states_modified<-c(top_states[top_states !="New York"],"Maryland")

# Plot with all states:
(Corona_Cases.US_state.summary.plot<-ggplot(Corona_Cases.US_state.summary,aes(x=Date,y=Total_confirmed_cases_perstate))+
  geom_point()+
  geom_point(data=filter(Corona_Cases.US_state.summary,Province.State %in% top_states),aes(col=Province.State))+
  scale_color_brewer(type = "qualitative",palette = "Paired")+
  default_theme+
  theme(axis.text.x = element_text(angle=45,hjust=1),legend.position = "bottom")+
  ggtitle("Total confirmed cases per state, top 12 colored")+
  scale_x_date(date_breaks="1 week",date_minor_breaks="1 day"))

##------------------------------------------
## Fit linear model to time vs total cases
##-----------------------------------------

# First, find the date at which each state's cases vs time becomes lienar (2nd derivative is about 0)
li<-ddply(Corona_Cases.US_state.summary,c("Province.State"),find_linear_index)

# Compute linear model for each state starting at the point at which data becomes linear
for(i in 1:nrow(li)){
  Province.State.i<-li[i,"Province.State"]
  date.i<-li[i,"V1"]
  data.i<-filter(Corona_Cases.US_state.summary,Province.State==Province.State.i & as.numeric(Date) >= date.i)
  model_results<-lm(data.i,formula = Total_confirmed_cases_perstate~Date)
  slope<-model_results$coefficients[2]
  intercept<-model_results$coefficients[1]
  li[li$Province.State==Province.State.i,"m"]<-slope
  li[li$Province.State==Province.State.i,"b"]<-intercept
  }

# Compute top state case load with fitted model

(Corona_Cases.US_state.lm.plot<-ggplot(filter(Corona_Cases.US_state.summary,Province.State %in% top_states_modified ))+
    geom_abline(data=filter(li,Province.State %in% top_states_modified),
                aes(slope = m,intercept = b,col=Province.State),lty=2)+
    geom_point(aes(x=Date,y=Total_confirmed_cases_perstate,col=Province.State))+
    scale_color_brewer(type = "qualitative",palette = "Paired")+
    default_theme+
    theme(axis.text.x = element_text(angle=45,hjust=1),legend.position = "bottom")+
    ggtitle("Total confirmed cases per state, top 12 colored")+
    scale_x_date(date_breaks="1 week",date_minor_breaks="1 day"))

##------------------------------------------
## Predict the number of total cases over the next week
##-----------------------------------------

predicted_days<-c(0,1,2,3,7)+as.numeric(as.Date("2020-04-20"))

predicted_days_df<-data.frame(matrix(ncol=3))
names(predicted_days_df)<-c("Province.State","days","Total_confirmed_cases_perstate")

# USe model parameters to estiamte case loads
for(state.i in top_states_modified){
  predicted_days_df<-rbind(predicted_days_df,
                           data.frame(Province.State=state.i,
                                      prediction_model(m = li[li$Province.State==state.i,"m"],
                                                       b =li[li$Province.State==state.i,"b"] ,
                                                       days =predicted_days )))
  }

predicted_days_df$Date<-as.Date(predicted_days_df$days,origin="1970-01-01")

kable(predicted_days_df,caption = "Predicted total cases over the next week for selected states")

##------------------------------------------
## Write plots
##-----------------------------------------

write_plot(Corona_Cases.US_state.summary.plot,wd = results_dir)
write_plot(Corona_Cases.US_state.lm.plot,wd = results_dir)

##------------------------------------------
## Write tables
##-----------------------------------------

write.csv(predicted_days_df,file = paste0(results_dir,"predicted_total_cases_days.csv"),quote = F,row.names = F)

Q2: What is the predicted number of cases?

What is the prediction of COVID-19 based on model thus far? Additional questions:

WHy did it take to day 40 to start a log linear trend? How long will it be till x number of cases? When will the plateu happen? Are any effects noticed with social distancing? Delays

##------------------------------------------
## Prediction and Prediction Accuracy
##------------------------------------------


today_num<-max(Corona_Cases.US$Days_since_100)
predicted_days<-today_num+c(1,2,3,7)

#mods = dlply(mydf, .(x3), lm, formula = y ~ x1 + x2)
#today:
Corona_Cases.US[Corona_Cases.US$Days_since_100==(today_num-1),]
Corona_Cases.US[Corona_Cases.US$Days_since_100==today_num,]
Corona_Cases.US$type<-"Historical"


#prediction_values<-prediction_model(m=slope,b=intercept,days = predicted_days)$Total_confirmed_cases

histoical_model<-data.frame(date=today_num,m=slope,b=intercept)
tmp<-data.frame(state=rep(c("A","B"),each=3),x=c(1,2,3,4,5,6))
tmp$y<-c(tmp[1:3,"x"]+5,tmp[4:6,"x"]*5+1)
ddply(tmp,c("state"))
lm(data =tmp,formula = y~x )

train_lm<-function(input_data,subset_coulmn,formula_input){
case_models <- dlply(input_data, subset_coulmn, lm, formula = formula_input)
case_models.parameters <- ldply(case_models, coef)
case_models.parameters<-rename(case_models.parameters,c("b"="(Intercept)","m"=subset_coulmn))
return(case_models.parameters)
}

train_lm(tmp,"state")

 dlply(input_data, subset_coulmn, lm,m=)
 
# model for previous y days
#historical_model_predictions<-data.frame(day_x=NULL,Days_since_100=NULL,Total_confirmed_cases=NULL,Total_confirmed_cases.log=NULL)
# for(i in c(1,2,3,4,5,6,7,8,9,10)){
#   #i<-1
# day_x<-today_num-i # 1, 2, 3, 4
# day_x_nextweek<-day_x+c(1,2,3)
# model_fit_x<-lm(data = filter(Corona_Cases.US.case100,Days_since_100 < day_x),formula = Total_confirmed_cases.log~Days_since_100)
# prediction_day_x_nextweek<-prediction_model(m = model_fit_x$coefficients[2],b = model_fit_x$coefficients[1],days = day_x_nextweek)
# prediction_day_x_nextweek$type<-"Predicted"
# acutal_day_x_nextweek<-filter(Corona_Cases.US,Days_since_100 %in% day_x_nextweek) %>% select(c(Days_since_100,Total_confirmed_cases,Total_confirmed_cases.log))
# acutal_day_x_nextweek$type<-"Historical"
# historical_model_predictions.i<-data.frame(day_x=day_x,rbind(acutal_day_x_nextweek,prediction_day_x_nextweek))
# historical_model_predictions<-rbind(historical_model_predictions.i,historical_model_predictions)
# }

#historical_model_predictions.withHx<-rbind.fill(historical_model_predictions,data.frame(Corona_Cases.US,type="Historical"))
#historical_model_predictions.withHx$Total_confirmed_cases.log2<-log(historical_model_predictions.withHx$Total_confirmed_cases,2)

(historical_model_predictions.plot<-ggplot(historical_model_predictions.withHx,aes(x=Days_since_100,y=Total_confirmed_cases.log,col=type))+
    geom_point(size=3)+
    default_theme+
    theme(legend.position = "bottom")+ 
      #geom_abline(slope = slope,intercept =intercept,lty=2)+
    #facet_wrap(~case_type,ncol=1)+
    scale_color_manual(values = c("Historical"="#377eb8","Predicted"="#e41a1c")))
write_plot(historical_model_predictions.plot,wd=results_dir)

Q3: What is the effect on social distancing, descreased mobility on case load?

Load data from Google which compoutes % change in user mobility relative to baseline for * Recreation
* Workplace
* Residence
* Park
* Grocery

Data from https://www.google.com/covid19/mobility/

# See pre-processing section for script on gathering mobility data

# UNDER DEVELOPMENT

mobility<-read.csv("/Users/stevensmith/Projects/MIT_COVID19/mobility.csv",header = T,stringsAsFactors = F)
#mobility$Retail_Recreation<-as.numeric(sub(mobility$Retail_Recreation,pattern = "%",replacement = ""))
#mobility$Workplace<-as.numeric(sub(mobility$Workplace,pattern = "%",replacement = ""))
#mobility$Residential<-as.numeric(sub(mobility$Residential,pattern = "%",replacement = ""))

##------------------------------------------
## Show relationship between mobility and caseload
##------------------------------------------
mobility$County<-gsub(mobility$County,pattern = " County",replacement = "")
Corona_Cases.US_state.mobility<-merge(Corona_Cases.US_state,plyr::rename(mobility,c("State"="Province.State","County"="City")))

#Corona_Cases.US_state.tmp<-merge(metadata,Corona_Cases.US_state.tmp)
# Needs to happen upsteam, see todos
#Corona_Cases.US_state.tmp$Total_confirmed_cases.perperson<-Corona_Cases.US_state.tmp$Total_confirmed_cases/as.numeric(Corona_Cases.US_state.tmp$Population)
mobility_measures<-c("Retail_Recreation","Grocery_Pharmacy","Parks","Transit","Workplace","Residential")

plot_data<-filter(Corona_Cases.US_state.mobility, Date.numeric==max(Corona_Cases.US_state$Date.numeric) ) %>% melt(measure.vars=mobility_measures) 
plot_data$value<-as.numeric(gsub(plot_data$value,pattern = "%",replacement = ""))
plot_data<-filter(plot_data,!is.na(value))

(mobility.plot<-ggplot(filter(plot_data,Province.State %in% c("Pennsylvania","Maryland","New Jersey","California","Delaware","Connecticut")),aes(y=Total_confirmed_cases.per100,x=value))+geom_point()+
  facet_grid(Province.State~variable,scales = "free")+
  xlab("Mobility change from baseline (%)")+
  ylab(paste0("Confirmed cases per 100 people(Today)"))+
  default_theme+
  ggtitle("Mobility change vs cases"))

(mobility.global.plot<-ggplot(plot_data,aes(y=Total_confirmed_cases.per100,x=value))+geom_point()+
  facet_wrap(~variable,scales = "free")+
  xlab("Mobility change from baseline (%)")+
  ylab(paste0("Confirmed cases (Today) per 100 people"))+
  default_theme+
  ggtitle("Mobility change vs cases"))

plot_data.permobility_summary<-ddply(plot_data,c("Province.State","variable"),summarise,cor=cor(y =Total_confirmed_cases.per100,x=value),median_change=median(x=value)) %>% arrange(-abs(cor))

kable(plot_data.permobility_summary,caption = "Ranked per-state mobility correlation with total confirmed cases")

Ranked per-state mobility correlation with total confirmed cases

Province.State	variable	cor	median_change
Alaska	Transit	-1.0000000	-63.0
Delaware	Retail_Recreation	1.0000000	-39.5
Delaware	Grocery_Pharmacy	1.0000000	-17.5
Delaware	Parks	-1.0000000	20.5
Delaware	Transit	1.0000000	-37.0
Delaware	Workplace	1.0000000	-37.0
Delaware	Residential	-1.0000000	14.0
Hawaii	Retail_Recreation	0.9931972	-56.0
Hawaii	Grocery_Pharmacy	0.9695437	-34.0
New Hampshire	Parks	0.9569819	-20.0
Connecticut	Grocery_Pharmacy	-0.9034799	-6.0
Maine	Transit	-0.9024399	-50.0
Alaska	Residential	0.8898278	13.0
South Dakota	Parks	0.8710801	-26.0
Utah	Residential	-0.8524968	12.0
Vermont	Parks	0.8519884	-35.5
Alaska	Grocery_Pharmacy	-0.8078182	-7.0
Hawaii	Residential	-0.7854909	19.0
Utah	Transit	-0.7794922	-18.0
Massachusetts	Workplace	-0.7711813	-39.0
Connecticut	Transit	-0.7510764	-50.0
Rhode Island	Workplace	-0.7500031	-39.5
Alaska	Workplace	-0.7296146	-34.0
Wyoming	Transit	-0.7287206	-17.0
Utah	Parks	-0.6985837	17.0
Wyoming	Parks	-0.6881358	-4.0
Hawaii	Parks	0.6813458	-72.0
Utah	Workplace	-0.6591668	-37.0
Vermont	Grocery_Pharmacy	-0.6585379	-25.0
New York	Workplace	-0.6458510	-34.5
Maine	Workplace	-0.6449558	-30.0
Rhode Island	Retail_Recreation	-0.6292286	-45.0
Montana	Workplace	-0.6239388	-40.5
Arizona	Grocery_Pharmacy	-0.6233132	-15.0
Hawaii	Transit	0.6188732	-89.0
Rhode Island	Residential	-0.6159667	18.5
New Jersey	Workplace	-0.6117849	-44.0
Nebraska	Workplace	0.6068569	-32.5
New Jersey	Parks	-0.5966373	-6.0
New York	Retail_Recreation	-0.5871506	-46.0
North Dakota	Retail_Recreation	-0.5596690	-42.0
Hawaii	Workplace	0.5396454	-46.0
Connecticut	Residential	0.5358815	14.0
Massachusetts	Retail_Recreation	-0.5293009	-44.0
New York	Parks	0.5268744	20.0
Connecticut	Retail_Recreation	-0.5165179	-45.0
New Jersey	Retail_Recreation	-0.5105718	-62.5
Maine	Parks	0.5013907	-31.0
Connecticut	Workplace	-0.4988833	-39.0
Arizona	Retail_Recreation	-0.4981312	-42.5
North Dakota	Parks	0.4973765	-34.0
Montana	Parks	-0.4913929	-58.0
New Jersey	Grocery_Pharmacy	-0.4858699	2.5
Nebraska	Residential	-0.4858344	14.0
New Mexico	Grocery_Pharmacy	-0.4793089	-11.0
Wyoming	Workplace	-0.4786493	-31.0
Iowa	Parks	-0.4777179	28.5
Rhode Island	Parks	0.4731974	52.0
Montana	Residential	0.4701424	14.0
New Mexico	Parks	0.4548477	-31.5
Illinois	Transit	-0.4517787	-31.0
New Mexico	Residential	0.4491212	13.5
Vermont	Residential	0.4379791	11.5
Idaho	Workplace	-0.4374387	-29.0
California	Transit	-0.4371164	-42.0
Kentucky	Parks	-0.4360494	28.5
California	Residential	0.4329382	14.0
Pennsylvania	Workplace	-0.4324064	-36.0
Wisconsin	Transit	-0.4318169	-23.5
Massachusetts	Grocery_Pharmacy	-0.4313347	-7.0
New Jersey	Transit	-0.4278509	-50.5
Kansas	Parks	0.4187728	72.0
South Carolina	Workplace	0.4180317	-30.0
Montana	Retail_Recreation	-0.4145442	-51.0
New Hampshire	Residential	-0.4124232	14.0
Arkansas	Parks	-0.4086800	-12.0
Idaho	Transit	-0.4084972	-30.0
Idaho	Grocery_Pharmacy	-0.4036292	-4.5
Alabama	Workplace	-0.3999512	-29.0
Maryland	Workplace	-0.3981259	-35.0
Arizona	Residential	0.3943607	13.0
Maryland	Grocery_Pharmacy	-0.3897732	-10.0
Montana	Transit	-0.3892824	-41.0
Alabama	Grocery_Pharmacy	-0.3814724	-2.0
New York	Transit	-0.3745865	-48.0
Arizona	Transit	0.3632915	-38.0
Florida	Residential	0.3592640	14.0
New Mexico	Retail_Recreation	-0.3574663	-42.5
Pennsylvania	Retail_Recreation	-0.3562258	-45.0
Wyoming	Grocery_Pharmacy	-0.3534260	-10.0
California	Parks	-0.3461554	-38.5
Alabama	Transit	-0.3436527	-36.5
Nevada	Transit	-0.3424973	-20.0
Michigan	Parks	0.3370423	30.0
Nebraska	Grocery_Pharmacy	0.3317520	-0.5
Pennsylvania	Parks	0.3292215	13.0
Montana	Grocery_Pharmacy	-0.3279939	-16.0
Alaska	Retail_Recreation	0.3260634	-39.0
West Virginia	Parks	0.3239758	-33.0
Minnesota	Transit	-0.3157689	-28.5
North Carolina	Grocery_Pharmacy	0.3100644	0.0
Idaho	Retail_Recreation	-0.3048434	-40.5
West Virginia	Grocery_Pharmacy	-0.3040290	-6.0
Maine	Retail_Recreation	-0.3010249	-42.0
Vermont	Retail_Recreation	0.2952799	-57.0
North Dakota	Workplace	0.2903163	-40.0
Colorado	Residential	0.2883815	14.0
Texas	Workplace	0.2807795	-32.0
California	Retail_Recreation	-0.2801884	-44.0
Utah	Retail_Recreation	-0.2800449	-40.0
Mississippi	Residential	0.2786968	13.0
Texas	Residential	-0.2786580	15.0
Arkansas	Retail_Recreation	-0.2768985	-30.0
Maryland	Retail_Recreation	-0.2758363	-39.0
Rhode Island	Transit	-0.2752347	-56.0
Kansas	Workplace	0.2739561	-32.5
California	Grocery_Pharmacy	-0.2723915	-11.5
Maryland	Residential	0.2682304	15.0
Oregon	Grocery_Pharmacy	0.2651842	-7.0
Florida	Parks	-0.2650683	-43.0
North Carolina	Workplace	0.2648515	-31.0
Nevada	Residential	0.2634265	17.0
California	Workplace	-0.2620151	-36.0
Nevada	Retail_Recreation	-0.2612579	-43.0
Virginia	Transit	-0.2600154	-33.0
Rhode Island	Grocery_Pharmacy	0.2565848	-7.5
Illinois	Workplace	-0.2525444	-31.0
Texas	Parks	0.2522483	-42.0
Tennessee	Workplace	-0.2513115	-31.0
Vermont	Workplace	-0.2476498	-43.0
Tennessee	Residential	0.2469408	11.5
Wisconsin	Parks	0.2459968	51.5
Georgia	Grocery_Pharmacy	-0.2451624	-10.0
South Carolina	Parks	-0.2442901	-23.0
Illinois	Parks	0.2393778	26.5
Pennsylvania	Grocery_Pharmacy	-0.2367309	-6.0
New York	Grocery_Pharmacy	-0.2318705	8.0
Arkansas	Residential	0.2316479	12.0
Washington	Workplace	-0.2199559	-38.0
Michigan	Workplace	-0.2171009	-40.0
Missouri	Residential	-0.2140702	13.0
North Carolina	Transit	0.2137332	-32.0
North Carolina	Residential	0.2119263	13.0
New Jersey	Residential	0.2104578	18.0
Missouri	Workplace	0.2011787	-28.5
Kansas	Grocery_Pharmacy	-0.1995442	-14.0
Iowa	Transit	0.1974885	-24.0
Oregon	Residential	0.1928806	10.5
South Dakota	Transit	-0.1921556	-40.0
Georgia	Workplace	-0.1891422	-33.5
Illinois	Residential	0.1870788	14.0
Mississippi	Grocery_Pharmacy	-0.1840467	-8.0
North Dakota	Grocery_Pharmacy	-0.1836997	-8.0
Idaho	Residential	-0.1806540	11.0
Colorado	Parks	-0.1791649	2.0
Wyoming	Retail_Recreation	-0.1771708	-39.0
Alabama	Parks	0.1718818	-1.0
Texas	Transit	0.1693244	-41.5
Georgia	Retail_Recreation	-0.1692230	-41.0
New Mexico	Transit	0.1679544	-38.5
Virginia	Grocery_Pharmacy	-0.1645141	-8.0
Ohio	Transit	0.1630571	-28.0
Wisconsin	Residential	-0.1629294	14.0
Georgia	Residential	-0.1618675	13.0
Virginia	Residential	0.1604840	14.0
Virginia	Parks	0.1596643	6.0
Oklahoma	Residential	0.1567643	15.0
South Carolina	Residential	-0.1538989	12.0
South Dakota	Retail_Recreation	-0.1506499	-38.5
Washington	Transit	-0.1504574	-33.5
Connecticut	Parks	0.1474430	43.0
Florida	Retail_Recreation	0.1465621	-43.0
Washington	Residential	0.1462162	13.0
Massachusetts	Transit	-0.1460918	-45.0
North Dakota	Transit	0.1433902	-48.0
Michigan	Retail_Recreation	-0.1432653	-53.0
New Hampshire	Retail_Recreation	-0.1429439	-41.0
Indiana	Retail_Recreation	0.1394149	-38.0
Indiana	Residential	0.1387390	12.0
Massachusetts	Parks	0.1379458	39.0
Maine	Residential	-0.1376179	11.0
Mississippi	Transit	-0.1373815	-38.5
North Carolina	Parks	-0.1315740	7.0
Florida	Workplace	-0.1314709	-33.0
Oregon	Retail_Recreation	0.1305017	-41.0
Pennsylvania	Transit	-0.1302178	-41.5
Alabama	Retail_Recreation	0.1295806	-39.0
South Dakota	Residential	0.1255465	15.0
Oklahoma	Parks	-0.1252441	-18.5
Minnesota	Parks	0.1245200	-9.0
Massachusetts	Residential	0.1227471	15.0
Ohio	Parks	-0.1206226	67.5
New Hampshire	Grocery_Pharmacy	-0.1185951	-6.0
Wisconsin	Workplace	-0.1183345	-31.0
Kansas	Transit	-0.1168930	-26.5
Maine	Grocery_Pharmacy	-0.1167170	-13.0
Washington	Grocery_Pharmacy	0.1163598	-7.0
Minnesota	Workplace	-0.1149393	-33.0
Oregon	Parks	0.1143547	16.5
Idaho	Parks	0.1110882	-22.0
Maryland	Transit	-0.1101439	-39.0
Mississippi	Retail_Recreation	-0.1099643	-40.0
Kentucky	Grocery_Pharmacy	0.1098535	4.0
Arkansas	Transit	0.1090221	-27.0
West Virginia	Residential	-0.1087555	11.0
Arkansas	Workplace	-0.1080883	-26.0
Ohio	Residential	0.1077552	14.0
Indiana	Parks	-0.1023437	29.0
New Hampshire	Transit	-0.1019342	-57.0
Arizona	Workplace	-0.1006545	-35.0
Nebraska	Retail_Recreation	0.0992677	-36.0
Minnesota	Retail_Recreation	0.0973131	-40.0
Mississippi	Workplace	-0.0956580	-33.0
Wisconsin	Grocery_Pharmacy	0.0954311	-1.0
Georgia	Parks	0.0951224	-6.0
Michigan	Grocery_Pharmacy	-0.0931909	-11.0
South Dakota	Workplace	0.0927870	-35.0
Texas	Grocery_Pharmacy	0.0923568	-14.0
Wyoming	Residential	0.0920486	12.5
Pennsylvania	Residential	0.0915801	15.0
New York	Residential	0.0913396	17.5
Nebraska	Parks	0.0911200	55.5
Oklahoma	Grocery_Pharmacy	-0.0895337	-1.0
Missouri	Transit	-0.0881870	-24.5
Indiana	Workplace	0.0871502	-34.0
Virginia	Workplace	-0.0869857	-31.5
Washington	Parks	0.0859680	-3.5
South Dakota	Grocery_Pharmacy	0.0843545	-9.0
Kentucky	Transit	0.0794671	-31.0
Michigan	Residential	0.0779151	15.0
Indiana	Grocery_Pharmacy	-0.0766021	-5.5
South Carolina	Transit	0.0760657	-45.0
Kentucky	Retail_Recreation	0.0747466	-29.0
Virginia	Retail_Recreation	-0.0742280	-35.0
Colorado	Transit	0.0719691	-36.0
Ohio	Grocery_Pharmacy	0.0666812	0.0
North Carolina	Retail_Recreation	0.0662733	-34.0
Nevada	Workplace	0.0655054	-40.0
Michigan	Transit	0.0651504	-46.0
Tennessee	Parks	-0.0645235	10.5
Ohio	Retail_Recreation	0.0622588	-36.0
West Virginia	Workplace	0.0612386	-33.0
Minnesota	Grocery_Pharmacy	0.0599051	-6.0
Oregon	Transit	0.0597446	-27.5
Nevada	Parks	0.0594394	-12.5
Oregon	Workplace	-0.0589402	-31.0
Nebraska	Transit	-0.0576652	-9.0
Iowa	Retail_Recreation	-0.0565372	-38.0
West Virginia	Retail_Recreation	-0.0526565	-38.5
Florida	Transit	-0.0516869	-49.0
Oklahoma	Workplace	0.0512421	-31.0
South Carolina	Retail_Recreation	-0.0496295	-35.0
South Carolina	Grocery_Pharmacy	0.0488516	1.0
Missouri	Parks	0.0484492	0.0
Colorado	Grocery_Pharmacy	-0.0482725	-17.0
Colorado	Retail_Recreation	-0.0480970	-44.0
Washington	Retail_Recreation	-0.0480542	-42.0
Illinois	Retail_Recreation	0.0475648	-40.0
Kentucky	Residential	0.0473650	12.0
Texas	Retail_Recreation	0.0470931	-40.0
New Hampshire	Workplace	0.0452092	-37.0
Kentucky	Workplace	-0.0449154	-36.0
Missouri	Grocery_Pharmacy	0.0441849	2.0
Missouri	Retail_Recreation	-0.0437999	-36.0
Iowa	Workplace	-0.0423004	-30.0
North Dakota	Residential	-0.0403637	17.0
Illinois	Grocery_Pharmacy	-0.0379889	2.0
Tennessee	Transit	-0.0376040	-32.0
Florida	Grocery_Pharmacy	0.0335606	-14.0
Arizona	Parks	-0.0322416	-44.5
Indiana	Transit	0.0291950	-29.0
Minnesota	Residential	-0.0288864	17.0
Vermont	Transit	0.0270377	-63.0
Tennessee	Grocery_Pharmacy	0.0262641	6.0
Ohio	Workplace	-0.0258992	-35.0
Oklahoma	Retail_Recreation	0.0230294	-31.0
Nevada	Grocery_Pharmacy	0.0189860	-12.5
West Virginia	Transit	-0.0170546	-45.0
Kansas	Residential	-0.0166925	13.0
Wisconsin	Retail_Recreation	0.0156966	-44.0
Mississippi	Parks	-0.0145460	-25.0
Kansas	Retail_Recreation	-0.0131605	-37.0
Utah	Grocery_Pharmacy	0.0130331	-4.0
Georgia	Transit	-0.0124298	-35.0
Maryland	Parks	-0.0105794	27.0
Iowa	Residential	-0.0081391	13.0
Colorado	Workplace	-0.0079635	-39.0
New Mexico	Workplace	0.0071133	-34.0
Oklahoma	Transit	0.0045994	-26.0
Tennessee	Retail_Recreation	-0.0036535	-30.0
Iowa	Grocery_Pharmacy	0.0025521	4.0
Arkansas	Grocery_Pharmacy	0.0024220	3.0
Alabama	Residential	-0.0009927	11.0
Alaska	Parks	NA	29.0
District of Columbia	Retail_Recreation	NA	-69.0
District of Columbia	Grocery_Pharmacy	NA	-28.0
District of Columbia	Parks	NA	-65.0
District of Columbia	Transit	NA	-69.0
District of Columbia	Workplace	NA	-48.0
District of Columbia	Residential	NA	17.0

# sanity check
ggplot(filter(plot_data,Province.State %in% c("Pennsylvania","Maryland","New Jersey","California","Delaware","Connecticut")),aes(x=Total_confirmed_cases.per100,fill=variable))+geom_histogram()+
  facet_grid(~Province.State)+
    default_theme+
  theme(legend.position = "bottom")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

write_plot(mobility.plot,wd = results_dir)

## [1] "/Users/stevensmith/Projects/coronavirus/results/mobility.plot.png"

write_plot(mobility.global.plot,wd = results_dir)

## [1] "/Users/stevensmith/Projects/coronavirus/results/mobility.global.plot.png"

(plot_data.permobility_summary.plot<-ggplot(plot_data.permobility_summary,aes(x=variable,y=median_change))+
  geom_jitter(size=2,width=.2)+
  #geom_jitter(data=plot_data.permobility_summary %>% arrange(-abs(median_change)) %>% head(n=15),aes(col=Province.State),size=2,width=.2)+
  default_theme+
  ggtitle("Per-Sate Median Change in Mobility")+
  xlab("Mobility Meaure")+
  ylab("Median Change from Baseline"))

write_plot(plot_data.permobility_summary.plot,wd = results_dir)

## [1] "/Users/stevensmith/Projects/coronavirus/results/plot_data.permobility_summary.plot.png"

DELIVERABLE MANIFEST

The following link to commited documents pushed to github. These are provided as a convienence, but note this is a manual process. The generation of reports, plots and tables is not coupled to the execution of this markdown. ## Report This report, html & pdf

Plots

github_root<-"https://github.com/sbs87/coronavirus/blob/master/"

plot_handle<-c("Corona_Cases.world.long.plot",
               "Corona_Cases.world.loglong.plot",
               "Corona_Cases.world.mortality.plot",
               "Corona_Cases.world.casecor.plot",
               "Corona_Cases.city.long.plot",
               "Corona_Cases.city.loglong.plot",
               "Corona_Cases.city.mortality.plot",
               "Corona_Cases.city.casecor.plot",
               "Corona_Cases.city.long.normalized.plot",
               "Corona_Cases.US_state.lm.plot",
               "Corona_Cases.US_state.summary.plot")


deliverable_manifest<-data.frame(
  name=c("World total & death cases, longitudinal",
         "World log total & death cases, longitudinal",
         "World mortality",
         "World total & death cases, correlation",
         "City total & death cases, longitudinal",
         "City log total & death cases, longitudinal",
         "City mortality",
         "City total & death cases, correlation",
         "City population normalized total & death cases, longitudinal",
         "State total cases (select) with linear model, longitudinal",
         "State total cases, longitudinal"),
  plot_handle=plot_handle,
  link=paste0(github_root,"results/",plot_handle,".png")
)


(tmp<-data.frame(row_out=apply(deliverable_manifest,MARGIN = 1,FUN = function(x) paste(x[1],x[2],x[3],sep=" | "))))

##                                                                                                                                                                                                        row_out
## 1                                           World total & death cases, longitudinal | Corona_Cases.world.long.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.long.plot.png
## 2                                 World log total & death cases, longitudinal | Corona_Cases.world.loglong.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.loglong.plot.png
## 3                                                         World mortality | Corona_Cases.world.mortality.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.mortality.plot.png
## 4                                      World total & death cases, correlation | Corona_Cases.world.casecor.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.casecor.plot.png
## 5                                              City total & death cases, longitudinal | Corona_Cases.city.long.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.long.plot.png
## 6                                    City log total & death cases, longitudinal | Corona_Cases.city.loglong.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.loglong.plot.png
## 7                                                            City mortality | Corona_Cases.city.mortality.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.mortality.plot.png
## 8                                         City total & death cases, correlation | Corona_Cases.city.casecor.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.casecor.plot.png
## 9  City population normalized total & death cases, longitudinal | Corona_Cases.city.long.normalized.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.long.normalized.plot.png
## 10                     State total cases (select) with linear model, longitudinal | Corona_Cases.US_state.lm.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.US_state.lm.plot.png
## 11                                      State total cases, longitudinal | Corona_Cases.US_state.summary.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.US_state.summary.plot.png

row_out<-apply(tmp, 2, paste, collapse="\t\n")

name	handle	link
World total & death cases, longitudinal	Corona_Cases.world.long.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.long.plot.png
World log total & death cases, longitudinal	Corona_Cases.world.loglong.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.loglong.plot.png
World mortality	Corona_Cases.world.mortality.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.mortality.plot.png
World total & death cases, correlation	Corona_Cases.world.casecor.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.casecor.plot.png
City total & death cases, longitudinal	Corona_Cases.city.long.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.long.plot.png
City log total & death cases, longitudinal	Corona_Cases.city.loglong.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.loglong.plot.png
City mortality	Corona_Cases.city.mortality.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.mortality.plot.png
City total & death cases, correlation	Corona_Cases.city.casecor.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.casecor.plot.png
City population normalized total & death cases, longitudinal	Corona_Cases.city.long.normalized.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.long.normalized.plot.png
State total cases (select) with linear model, longitudinal	Corona_Cases.US_state.lm.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.US_state.lm.plot.png
State total cases, longitudinal	Corona_Cases.US_state.summary.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.US_state.summary.plot.png

Tables

CONCLUSION

Overall, the trends of COVID-19 cases is no longer in log-linear phase for world or U.S. (but some regions like MD are still in the log-linear phase). Mortality rate (deaths/confirmed RNA-based cases) is >1%, with a range depending on region. Mobility is not a strong indicator of caseload (U.S. data).

See table below for detailed breakdown.

Question	Answer
What is the effect on social distancing, descreased mobility on case load?	There is not a strong apparent effect on decreased mobility (work, grocery, retail) or increased mobility (at residence, parks) on number of confirmed cases, either as a country (U.S.) or state level. California appears to have one of the best correlations, but this is a mixed bag
What is the trend in cases, mortality across geopgraphical regions?	The confirmed total casees and mortality is overall log-linear for most countries, with a trailing off beginning for most (inlcuding U.S.). On the state level, NY, NJ, PA starting to trail off; MD is still in log-linear phase. Mortality and case load are highly correlated for NY, NJ, PA, MD. The mortality rate flucutates for a given region, but is about 3% overall.

END

End: ##—— Sun May 24 10:47:15 2020 ——##

Cheatsheet: http://rmarkdown.rstudio.com>

Sandbox

# Geographical heatmap!
install.packages("maps")
library(maps)
library
mi_counties <- map_data("county", "pennsylvania") %>% 
  select(lon = long, lat, group, id = subregion)
head(mi_counties)

ggplot(mi_counties, aes(lon, lat)) + 
  geom_point(size = .25, show.legend = FALSE) +
  coord_quickmap()
mi_counties$cases<-1:2226
name_overlaps(metadata,Corona_Cases.US_state)

tmp<-merge(Corona_Cases.US_state,metadata)
ggplot(filter(tmp,Province.State=="Pennsylvania"), aes(Long, Lat, group = as.factor(City))) +
  geom_polygon(aes(fill = Total_confirmed_cases), colour = "grey50") + 
  coord_quickmap()


ggplot(Corona_Cases.US_state, aes(Long, Lat))+
  geom_polygon(aes(fill = Total_confirmed_cases ), color = "white")+
  scale_fill_viridis_c(option = "C")
dev.off()


require(maps)
require(viridis)

world_map <- map_data("world")
ggplot(world_map, aes(x = long, y = lat, group = group)) +
  geom_polygon(fill="lightgray", colour = "white")

head(world_map)
head(Corona_Cases.US_state)
unique(select(world_map,c("region","group"))) %>% filter()

some.eu.countries <- c(
  "US"
)
# Retrievethe map data
some.eu.maps <- map_data("world", region = some.eu.countries)

# Compute the centroid as the mean longitude and lattitude
# Used as label coordinate for country's names
region.lab.data <- some.eu.maps %>%
  group_by(region) %>%
  summarise(long = mean(long), lat = mean(lat))

unique(filter(some.eu.maps,subregion %in% Corona_Cases.US_state$Province.State) %>% select(subregion))
unique(Corona_Cases.US_state$Total_confirmed_cases.log)
ggplot(filter(Corona_Cases.US_state,Date=="2020-04-17") aes(x = Long, y = Lat)) +
  geom_polygon(aes( fill = Total_confirmed_cases.log))+
  #geom_text(aes(label = region), data = region.lab.data,  size = 3, hjust = 0.5)+
  #scale_fill_viridis_d()+
  #theme_void()+
  theme(legend.position = "none")
library("sf")
library("rnaturalearth")
library("rnaturalearthdata")

world <- ne_countries(scale = "medium", returnclass = "sf")
class(world)
ggplot(data = world) +
    geom_sf()

counties <- st_as_sf(map("county", plot = FALSE, fill = TRUE))
counties <- subset(counties, grepl("florida", counties$ID))
counties$area <- as.numeric(st_area(counties))
#install.packages("lwgeom")
class(counties)
head(counties)
ggplot(data = world) +
    geom_sf(data=Corona_Cases.US_state) +
    #geom_sf(data = counties, aes(fill = area)) +
  geom_sf(data = counties, aes(fill = area)) +
   # scale_fill_viridis_c(trans = "sqrt", alpha = .4) +
    coord_sf(xlim = c(-88, -78), ylim = c(24.5, 33), expand = FALSE)


head(counties)
tmp<-unique(select(filter(Corona_Cases.US_state,Date=="2020-04-17"),c(Lat,Long,Total_confirmed_cases.per100)))
st_as_sf(map("county", plot = FALSE, fill = TRUE))

join::inner_join.sf(Corona_Cases.US_state, counties)

library(sf)
library(sp)

nc <- st_read(system.file("shape/nc.shp", package="sf"))
class(nc)


spdf <- SpatialPointsDataFrame(coords = select(Corona_Cases.US_state,c("Lat","Long")), data = Corona_Cases.US_state,
                               proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))

head(spdf)
class(spdf)
st_cast(spdf)

filter(Corona_Cases.US_state.summary,Date=="2020-04-20" & Province.State %in% top_states_modified)

https://stevenbsmith.net